Classifying and reducing errors in density functional calculations 
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The energy error of any DFT calculation can be decomposed into a contribution due to the 
approximate functional and that due to the approximate density. Typically, the functional error 
dominates, but in many interesting situations, the density-driven error dominates. Examples include 
electron affinities, dissociation of molecules into charged fragments, transition state barriers, and 
ions and radicals in solution. In these abnormal cases, the error can be greatly reduced by using 
a more accurate density. A small gap between occupied and unoccupied orbitals may indicate a 
substantial density-driven error. 



Each year, at least 10,000 calculations of the elec- 
tronic structure of atoms, molecules, and solids, using 
Kohn-Sham (KS) density functional theory (DFT) pQ, 
are reported[5j. In all such calculations, a small but vital 
part of the total energy, called the exchange-correlation 
(XC) energy, is approximated. In any DFT calculation, 
one can consider the error as arising from two distinct 
sources. The Euler equation minimizing the total energy 
gives rise to a self-consistent equation for the density in- 
volving the functional derivative of the energy^. In KS- 
DFT, these are the celebrated KS equations PQ and in- 
clude the XC potential. Insertion of the solution density 
into the original energy functional yields the approximate 
energy. We write the error in an energy as the sum of 
two terms: 

AE — AEp + AEp,. (1) 

The functional error, AEp, is the contribution due to 
approximating the energy functional, while AEp, is the 
density-driven error due to the error in the self-consistent 
density. Under typical circumstances of KS calculations, 
self-consistent densities are extremely accurate, and the 
error in such normal calculations is dominated by AEp. 
The functional error of well-known approximations for 
given properties of given systems has become known from 
experience^, and can often be rationalized [5 ]. More 
sophisticated approximations are often tested using the 
densities (and orbitals) from less-sophisticated ones [6], a 
practice that assumes the calculation is normal. 

But in a small fraction (but still a very large number) 
of calculations, AEp, dominates over AEp, These errors 
can be greatly reduced by using more accurate densities. 
More importantly, the large energy error is typically mis- 
classified as being due to a limitation of the functional 
approximation, rather than the abnormality of the calcu- 
lation. The infamous self-interaction error (SIE) [7] made 
by standard DFT approximations for small atomic anions 
is so large that the ground-state is not even boundOE]. 
In fact, these are abnormal calculations and electron 
affinities with standard approximations evaluated on bet- 
ter densities are extremely accurate |101 HTj . 
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FIG. 1. Errors in ground-state energies of two-electron ions 
as a function of nuclear charge: PBE energies evaluated on 
exact [T2"lll3| (solid) and Hartree-Fock (HF, dotted) densities. 

Fig. [T] shows the paradigm of a catastrophic density- 
driven error, the two-electron ion as a function of Z, the 
nuclear charge. As has long been known [L4"j. as Z is 
reduced from 2 (He) to 1 (H~), a fraction of an electron 
unbinds (about 0.3) in a standard DFT calculation|15j. 
leading to an increasing error with decreasing Z (solid 
black line). But the solid colored lines show that, around 
Z c rj 1.23, where the orbital energy vanishes and the 
system begins to ionize, it is the density-driven error that 
is growing, and leads to the qualitative change in AE. 
The functional error is almost zero for H _ and is far less 
than for He. A DFT calculation with an accurate two- 
electron density produces a smaller error for the electron 
affinity of H than for the ionization energy of Hc[10 . 

In this letter, we show that (i) abnormality is a general 
phenomenon that can affect broad classes of DFT calcu- 
lations, from two-electron ions to radicals in solution, (ii) 
more accurate densities qualitatively improve such calcu- 
lations, (iii) a small (or zero) HOMO-LUMO gap in the 
approximate DFT calculation is an indicator of a density- 
driven error, and (iv) HF densities are sufficient to cure 
many density-driven errors. 
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We begin with a more detailed definition of Eq. (1). 
Let E be the exact energy of a system and n(r) its elec- 
tronic density. Let E and n(r) be their analogs in the 
approximate calculation. Then 

AEp = E[n] — E[n], AE D = E[h] - E[n\. (2) 

Since only the XC contribution is approximated in KS- 
DFT, AE F — E xo [n] — E xc [n\. We apply our definition 
to both total energies and energy differences. 
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FIG. 2. Exact and PBE radial densities for H~ and He. 

Most DFT calculations use either a generalized gradi- 
ent approximation (GGA) [T5] or a hybrid functional [151 
117] . The functional error is transmitted to the density via 
the XC potential, v xc (r) = SE XC /Sn(r). For these stan- 
dard approximations, which contain SIE, a well-known 
failing is that v xc (r) decays too rapidly with r, so that 
the overall KS potential is too shallow[T5]. For normal 
calculations, the KS eigenvalues are insufficiently deep, 
by several eV. However, a constant shift in the potential 
has no effect on n(r), which is very accurate, as seen in 
the inset of Fig. [2] This is quantified by the magnitude 
of AE D in Fig. [I] for Z = 2. Thus, normally, the SIE in 
the potential produces little error in the total energy. 

In an abnormal calculation, the system is peculiarly 
sensitive to the error in v xc (r), so that the n(r) differs 
significantly from n(r), enhancing AEe>. To make Fig. [T 
we carefully interpolated accurate QMC energies [121 ITS" 
from Z^ 1 = 1 to 0, and repeated this procedure apply- 
ing PBE to the exact densities to find AEp(Z). We also 
used Turbomole 19J to solve for self-consistent PBE en- 
ergies and the eigenvalue. For Z < Z c , the eigenvalue 
is pinned to 0, and an increasing fraction of an electron 
escapes. To achieve self consistency, we decrease the oc- 
cupation of the orbital from 2 until we find an occupation 
at which the total energy converges and the eigenvalue 
vanishes. The large error in density is very visible in 
Fig. [2] where the PBE density integrates to only 1.7 
electrons, so that all energy components, especially the 



Hartrcc energy, are quite inaccurate. Because of the vari- 
ational principle, the total energy remains more accurate 
than any of the individual components. From Fig. [T] one 
sees that PBE is (somewhat accidentally) almost exact 
for H~ when evaluated on the exact density. This illus- 
trates the basic difference between normal and abnormal 
calculations: Normally, improving the density has only a 
negligible effect on performance; in an abnormal calcula- 
tion, it can qualitatively improve results. 

Our method for classifying DFT errors is general. 
Declaring a calculation abnormal depends on both the 
energy being calculated (total, ionization, bond, etc.) 
and the approximation being used. The error in any 
approximation can be studied in this way. We focus here 
on the SIE of standard approximate functionals because 
of its ubiquity, but one can apply the same reasoning to, 
e.g., the correlation error in optimized effective potential 
(OEP) exact exchange calculations [20] , or the error in 
the KS kinetic energy in orbital- free approximations [2Tj . 
Much recent research is focussed on localization errors 
of approximations [22]. The classic examples of stretched 
H2 and H^" are normal, because self-consistent densities 
(restricted in the case of H2) are not so different from 
exact densities. More complex situations should be re- 
examined. 

We do not always need an extremely accurate density 
to eliminate density-driven errors, merely a density that 
is substantially more accurate than the abnormal density 
of the DFT calculation. For atomic and small molec- 
ular anions, a simple HF calculation is often sufficient 
to largely eliminate the density-driven error of standard 
DFT approximations. DFT energies evaluated on HF 
densities, called HF-DFT[TU], yield more accurate elec- 
tron affinities than DFT ionization potentials(Hj. 

An indicator of a likely density-driven error is given by 
the HOMO-LUMO gap in the approximate calculation. 
Apply linear response theory to the KS system: 

Sn(r) = J d 3 / Xs (r,r')^ s (r') (3) 

where 5n(r) is the change in density induced by 6v s (r), 

Xs (r, r ) = J2(fi " fj) lQ (4) 

is the static density-density KS response function, and 
4>i(r), £i, and fi are the KS orbitals, eigenvalues, and oc- 
cupation factors, respectively [3]. Since the occupation 
factors determine that i only includes occupied orbitals 
while j sums over only unoccupied orbitals, the small- 
est denominator (and hence the largest contribution) is 
from the transition from the highest occupied KS orbital 
to the lowest unoccupied one, called the HOMO-LUMO 
transition. Normally, the difference between the exact 
and approximate KS potentials is small, ignoring any 
constant shift. If the approximate HOMO-LUMO gap 
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(Ae g ) is not unusually small, this error leads to a small 
error in density. But if Ae g is small, even a small er- 
ror in the KS potential can produce a large change in 
the density, and self consistency only increases this ef- 
fect. Thus small Ae g suggest large density errors, and 
we include the gap in Fig. [I] For two-electron ions, 
the PBE LUMO is unbound, so that Ae g = \enoMo\- 
At Z c this vanishes. For atomic anions in general, the 
HOMO of the PBE potential on the exact density has a 
positive energy, i.e., a resonance [TO]- Finite-basis atom- 
centered calculations turn this resonance into an eigen- 
state, making small density-driven errors, and produce 
accurate electron affinities! 23 . 
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FIG. 3. Energy of as a function of separation in several 
calculations, and the PBE HOMO-LUMO gap. 

However, a small gap does not always imply a density- 
driven error. Consider the classic example of a severe 
SIE, namely stretched Hj with a standard functional |22j. 
as shown in Fig. [3j The small gap suggests abnormality, 
but when a HF-DFT calculation is performed, the error 
barely changes. Thus, this is a normal calculation, these 
errors are functional-driven, not density-driven, and HF- 
DFT is not recommended. 

Our next abnormality is well-known |24| . DFT calcula- 
tions of molecular dissociation energies (Eb) are usefully 
accurate with GGA's, and more so with hybrid func- 
tionals. These errors are often about 0.1 eV /bond [2"5]. 
found by subtracting the calculated molecular energy at 
its minimum from the sum of calculated atomic energies. 
This is because, if one simply increases the bond lengths 
to very large values, the fragments fail to dissociate into 
neutral atoms. The prototypical case is NaCl, which dis- 
sociates into Na 4 and CI -0 4 in a PBE calculation [2~4"j. 
The large error in density for the stretched bond yields 
AEf, « 1 eV, as shown as the difference between PBE and 
HF-PBE in Fig. g] In this case, the HF density spon- 
taneously suddenly switches to neutral atoms at about 
5.6A, but is correct in the dissociation limit. The com- 
mon practice of using isolated atomic calculations is in- 



consistent, but removes the density-driven error, because 
isolated atoms are normal. Incorrect dissociation occurs 
whenever the approximate HOMO of one is below the 
LUMO of the other [M], which guarantees a vanishing 
Ae g when the bond is greatly stretched. The exact v xc (r) 
contains a step between the atoms that ensures correct 
dissociation, but which is missed by semilocal approxi- 
mations. 
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FIG. 4. Energy of NaCl as a function of Na-Cl distance in 
several calculations, and the PBE HOMO-LUMO gap. 

When modern functionals were first being adopted for 
molecular calculations, they were sometimes evaluated 
on HF densities [2"6"ll2"T] . so as to compare only functional 
errrors. More recently. Janesko and Scuseria[55] showed 
this led to substantial improvement in transition-state 
barriers. The prototype of such barriers is the symmet- 
ric H— H2 transition state, which is improved by almost 
a factor of 2 by using HF-PBE instead of PBE. Here 
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is not quite as small (2.5 eV) as in other cases, 



but the improvement upon using the HF densities is still 
substantial. High-level ab initio calculations yield an en- 
ergy barrier of 0.43 eV[28] , where PBE gives value of 0.16 
eV and HF-PBE gives 0.25 eV. Analysis of a collection of 
barriers in Table 1 of Ref. I2U1 shows that, in cases where 
the HF-DFT barrier differs from the self-consistent bar- 
rier by more than, for instance 25%, the mean absolute 
error is more than three times smaller than DFT, with 
the sole exception of LN2H2 hydrogen transfer toward 
reaction barrier, where HF density is spin-contaminated 
(just as in the molecule CNjllj). 

Finally we report new applications where we drive out 
the density-driven error. We investigate potential en- 
ergy surfaces (PES) of odd-electron radical complexes 
like OH-CU and OHH 2 0. These are important in ra- 
diation, atmospheric, and environmental chemistry, as 
well as in cell biology [30H33] , For example, the behav- 
ior of anions in droplets is critical to our understanding 
aerosols in atmospheric chemistry 34J . Accepted wisdom 
claims that anions near an air-water interface, being less 



screened, have lower concentrations [35] ■ But recent clas- 
sical molecular dynamics (MD) simulations have shown 
the opposite[35tl36]. This controversy invites an ab initio 
MD approach, to either reinforce or debunk the classical 
simulations. 

However, DFT approximations have problems with 
these systems [3TH4*0] . Several such PES studies show 
two minima in the ground-state PES. One is a hydrogen- 
bonding structure while the other is a 2-center, 3-electron 
interacting hemi-bonding structure|41|. High-level quan- 
tum chemical calculations |41j and self- interaction cor- 
rected DFT results gQl [42] reveal that the true PES has 
only one minimum, the hydrogen-bonding structure. The 
hemibonding structure is relatively overstabilized in ap- 
proximate DFT because three electrons are incorrectly 
delocalized over two atoms. 
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FIG. 5. Potential energy surface scans for the HO-C1 - com- 
plex along CI - O - H angle of 0° (hydrogen-bonding structure, 
black lines) and 60° (hemi-bonding structure, red lines) using 
various methods; R is the Cl-O separation. 

In Fig. [5j we show plots of the PES of HO -CI - complex 
using different methods. The O-H bond length was fixed 
at 1 A. The binding energy is 
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FIG. 6. Binding energies in eV of HOH2O calculated with 
various methods for (a) hydrogen-bonding structure and (b) 
hemi-bonding structure. 

PBE greatly overstabilizes the hemi-bonding region, re- 
sulting in a strong, unphysical hemi-bonding minimum. 
CCSD(T) and HF-PBE, on the other hand, give an in- 
tuitively appealing hydrogen-bonding geometry as the 
global minimum. But Ae g of PBE is less than I eV, 
and HF-PBE fixes the problem. 

To summarize, we have shown how all approximate 
DFT calculations can be analyzed, to check for density- 
driven errors. In the case of SIEs of commonly-used func- 
tionals, use of HF (or better) densities can often greatly 
reduce such errors. This would allow distinctions be- 
tween true SIE (i.e., large energy errors even on exact 
densities) versus density-driven errors, which are system- 
and property-dependent. 

Calculation details - For evaluation of DFT and HF 
self-consistent densities in Fig 1 and 2, aug-cc-pV6Z ba- 
sis set [33] was used. From Fig 3-6, all calcuations were 
performed with aug-cc-pVTZ[44-46] basis with Turbo- 
mole. 
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HO-Cl 



(R, 9) - {E HQ . + E, 



cr 



(5) 



where -Ehoci- (-^i ^) * s the energy on a given geometry 
with Cl-0 distance R, and CI - O - H angle 9, Euo- is 
the energy of the OH radical, and E CI - is the energy 
of CI anion. The difference between the energy min- 
ima of hydrogen- and hemi-bonding structures in PBE 
is less than 0.01 eV. A small Ae g of PBE in the hy- 
drogen bonding structure (~ 0.32 eV) suggests a large 
density-driven error. We find HF-PBE follows the same 
trends and produces the same mimina as CCSD(T), al- 
though the binding energies themselves have errors of up 
to 0.09 eV. We also scanned the PES of the HOH 2 com- 
plex. The minimum geometry of hydrogen- and hemi- 
bonding structures is shown in Fig. [6] Self-consistent 
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